home *** CD-ROM | disk | FTP | other *** search
/ APDL Eductation Resources / APDL Eductation Resources.iso / programs / electronic / rlab / TestMatrix / randsvd_r < prev    next >
Encoding:
Text File  |  1994-12-20  |  3.9 KB  |  145 lines
  1. //-------------------------------------------------------------------//
  2.  
  3. // Synopsis:    Random matrix with pre-assigned singular values.
  4.  
  5. // Syntax:    R = randsvd (N, KAPPA, MODE, KL, KU)
  6.  
  7. // Description:
  8.  
  9. //    R is a (banded) random matrix of order N with COND(A) = KAPPA
  10. //    and singular values from the distribution MODE.
  11.  
  12. //      N may be a 2-vector, in which case the matrix is N(1)-by-N(2).
  13. //      Available types:
  14. //             MODE = 1:   one large singular value,
  15. //             MODE = 2:   one small singular value,
  16. //             MODE = 3:   geometrically distributed singular values,
  17. //             MODE = 4:   arithmetically distributed singular values,
  18. //             MODE = 5:   random singular values with unif. dist. logarithm.
  19.  
  20. //      If omitted, MODE defaults to 3, and KAPPA defaults to SQRT(1/EPS).
  21. //      If MODE < 0 then the effect is as for ABS(MODE) except that in the
  22. //      original matrix of singular values the order of the diagonal entries
  23. //      is reversed: small to large instead of large to small.
  24. //      KL and KU are the lower and upper bandwidths respectively; if they
  25. //      are omitted a full matrix is produced.
  26. //      If only KL is present, KU defaults to KL.
  27. //      Special case: if KAPPA < 0 then a random full symmetric positive
  28. //                    definite matrix is produced with COND(A) = -KAPPA and
  29. //                    eigenvalues distributed according to MODE.
  30. //                    KL and KU, if present, are ignored.
  31.  
  32. //      This routine is similar to the more comprehensive Fortran routine xLATMS
  33. //      in the following reference:
  34. //      J.W. Demmel and A. McKenney, A test matrix generation suite,
  35. //      LAPACK Working Note #9, Courant Institute of Mathematical Sciences,
  36. //      New York, 1989.
  37.  
  38. //    This file is a translation of randsvd.m from version 2.0 of
  39. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  40. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  41.  
  42. // Dependencies
  43.    require bandred qmult
  44.  
  45. //-------------------------------------------------------------------//
  46.  
  47. randsvd = function (n, kappa, mode, kl, ku)
  48. {
  49.   local (n, kappa, mode, kl, ku, factor)
  50.   global (eps)
  51.  
  52.   if (!exist (kappa)) { kappa = sqrt(1/eps); }
  53.   if (!exist (mode))  { mode = 3; }
  54.   if (!exist (kl))    { kl = n-1; }    // Full matrix.
  55.   if (!exist (ku))    { ku = kl; }     // Same upper and lower bandwidths.
  56.  
  57.   if (abs(kappa) < 1) {
  58.     error("Condition number must be at least 1!");
  59.   }
  60.  
  61.   posdef = 0;
  62.   if (kappa < 0)            // Special case.
  63.   {
  64.     posdef = 1;
  65.     kappa = -kappa;
  66.   }
  67.  
  68.   p = min(n);
  69.   m = n[1];        // Parameter n specifies dimension: m-by-n.
  70.   n = n[max(size(n))];
  71.  
  72.   if (p == 1)        // Handle case where A is a vector.
  73.   {
  74.     rand("normal", 0, 1);
  75.     A = rand(m, n);
  76.     A = A/norm(A, "2");
  77.     return A;
  78.   }
  79.  
  80.   j = abs(mode);
  81.  
  82.   // Set up vector sigma of singular values.
  83.   if (j == 3)
  84.   {
  85.     factor = kappa^(-1/(p-1));
  86.     sigma = factor.^[0:p-1];
  87.  
  88.   else if (j == 4) {
  89.     sigma = ones(p,1) - (0:p-1)'/(p-1)*(1-1/kappa);
  90.  
  91.   else if (j == 5) {    // In this case cond(A) <= kappa.
  92.     rand("uniform", 0, 1)
  93.     sigma = exp( -rand(p,1)*log(kappa) );
  94.  
  95.   else if (j == 2) {
  96.     sigma = ones(p,1);
  97.     sigma[p] = 1/kappa;
  98.  
  99.   else if (j == 1) {
  100.     sigma = ones(p,1)./kappa;
  101.     sigma[1] = 1;
  102.  
  103.   }}}}}
  104.  
  105.  
  106.   // Convert to diagonal matrix of singular values.
  107.   if (mode < 0) {
  108.     sigma = sigma[p:1:-1];
  109.   }
  110.  
  111.   sigma = diag(sigma);
  112.  
  113.   if (posdef)        // Handle special case.
  114.   {
  115.     Q = qmult(p);
  116.     A = Q'*sigma*Q;
  117.     A = (A + A')/2;    // Ensure matrix is symmetric.
  118.     return A;
  119.   }
  120.  
  121.   if (m != n) 
  122.   {
  123.     sigma[m; n] = 0;    // Expand to m-by-n diagonal matrix.
  124.   }
  125.  
  126.   if (kl == 0 && ku == 0)    // Diagonal matrix requested - nothing more to do.
  127.   {
  128.     A = sigma;
  129.     return A;
  130.   }
  131.  
  132.   // A = U*sigma*V, where U, V are random orthogonal matrices from the
  133.   // Haar distribution.
  134.   A = qmult(sigma');
  135.   A = qmult(A');
  136.  
  137.   if (kl < n-1 || ku < n-1)    // Bandwidth reduction.
  138.   {
  139.    A = bandred(A, kl, ku);
  140.   }
  141.  
  142.   rand("default");
  143.   return A;
  144. };
  145.